#ifndef cathlibcpp_complex_H
#define cathlibcpp_complex_H

// File:       complex.h
// Author:     (c) Miles Sabin, 1996
// Purpose:    approximation to ANSI C++ complex template


#ifndef included_math_H
#define included_math_H
#include <math.h>
#endif

#ifndef cathlibcpp_bool_H
#include "bool.h"
#endif

#ifndef cathlibcpp_config_H
#include "config.h"
#endif

#ifndef cathlibcpp_iosfwd_H
#include "iosfwd.h"              // for istream, ostream
#endif


template<class T>
class complex
{
  public:

    // constructors
    inline complex();
    inline complex(T const& re);
    inline complex(T const& re, T const& im);
    inline complex(complex<float> const& rhs);
    inline complex(complex<double> const& rhs);

    // accessors
    inline T real() const;
    inline T imag() const;

    // mutators
    inline complex<T>& operator=(T const& rhs);
    inline complex<T>& operator=(complex<float> const& rhs);
    inline complex<T>& operator=(complex<double> const& rhs);

    inline complex<T>& operator+=(T const& rhs);
    inline complex<T>& operator+=(complex<float> const& rhs);
    inline complex<T>& operator+=(complex<double> const& rhs);

    inline complex<T>& operator-=(T const& rhs);
    inline complex<T>& operator-=(complex<float> const& rhs);
    inline complex<T>& operator-=(complex<double> const& rhs);

    inline complex<T>& operator*=(T const& rhs);
    complex<T>& operator*=(complex<float> const& rhs);
    complex<T>& operator*=(complex<double> const& rhs);

    inline complex<T>& operator/=(T const& rhs);
    complex<T>& operator/=(complex<float> const& rhs);
    complex<T>& operator/=(complex<double> const& rhs);

  private:

    T re_;
    T im_;
};


template<class T>
inline complex<T> operator+(complex<T> const& x, complex<T> const& y);

template<class T>
inline complex<T> operator+(complex<T> const& x, T const& y);

template<class T>
inline complex<T> operator+(T const& x, complex<T> const& y);

template<class T>
inline complex<T> operator-(complex<T> const& x, complex<T> const& y);

template<class T>
inline complex<T> operator-(complex<T> const& x, T const& y);

template<class T>
inline complex<T> operator-(T const& x, complex<T> const& y);

template<class T>
inline complex<T> operator*(complex<T> const& x, complex<T> const& y);

template<class T>
inline complex<T> operator*(complex<T> const& x, T const& y);

template<class T>
inline complex<T> operator*(T const& x, complex<T> const& y);

template<class T>
inline complex<T> operator/(complex<T> const& x, complex<T> const& y);

template<class T>
inline complex<T> operator/(complex<T> const& x, T const& y);

template<class T>
complex<T> operator/(T const& x, complex<T> const& y);

template<class T>
inline complex<T> operator+(complex<T> const& x);

template<class T>
inline complex<T> operator-(complex<T> const& x);

template<class T>
inline bool operator==(complex<T> const& x, complex<T> const& y);

template<class T>
inline bool operator==(complex<T> const& x, T const& y);

template<class T>
inline bool operator==(T const& x, complex<T> const& y);

template<class T>
inline bool operator!=(complex<T> const& x, complex<T> const& x);

template<class T>
inline bool operator!=(complex<T> const& x, T const& y);

template<class T>
inline bool operator!=(T const& x, complex<T> const& y);

template<class T>
basic_istream_char& operator>>(basic_istream_char& is, complex<T>& x);

template<class T>
basic_ostream_char& operator<<(basic_ostream_char& os, complex<T> const& x);

template<class T>
inline T real(complex<T> const& x);

template<class T>
inline T imag(complex<T> const& x);

template<class T>
inline T abs(complex<T> const& x);

template<class T>
inline T arg(complex<T> const& x);

template<class T>
inline T norm(complex<T> const& x);

template<class T>
inline complex<T> conj(complex<T> const& x);

template<class T>
inline complex<T> polar(T const& mag, T const& theta);

template<class T>
inline complex<T> cos(complex<T> const& x);

template<class T>
inline complex<T> cosh(complex<T> const& x);

template<class T>
inline complex<T> exp(complex<T> const& x);

template<class T>
inline complex<T> log(complex<T> const& x);

template<class T>
inline complex<T> log10(complex<T> const& x);

template<class T>
complex<T> pow(complex<T> const& x, int y);

template<class T>
inline complex<T> pow(complex<T> const& x, T const& y);

template<class T>
inline complex<T> pow(complex<T> const& x, complex<T> const& y);

template<class T>
inline complex<T> pow(T const& x, complex<T> const& y);

template<class T>
inline complex<T> sin(complex<T> const& x);

template<class T>
inline complex<T> sinh(complex<T> const& x);

template<class T>
complex<T> sqrt(complex<T> const& x);

template<class T>
inline complex<T> tan(complex<T> const& x);

template<class T>
inline complex<T> tanh(complex<T> const& x);


// Implementation of complex<T>

template<class T>
inline complex<T>::complex()
  {}

template<class T>
inline complex<T>::complex(T const& re)
  : re_(re),
    im_(0)
  {}

template<class T>
inline complex<T>::complex(T const& re, T const& im)
  : re_(re),
    im_(im)
  {}

template<class T>
inline complex<T>::complex(complex<float> const& rhs)
  : re_(rhs.real()),
    im_(rhs.imag())
  {}

template<class T>
inline complex<T>::complex(complex<double> const& rhs)
  : re_(rhs.real()),
    im_(rhs.imag())
  {}

template<class T>
inline T complex<T>::real() const
  { return re_; }

template<class T>
inline T complex<T>::imag() const
  { return im_; }

template<class T>
inline complex<T>& complex<T>::operator=(T const& rhs)
  {
    re_ = rhs;
    im_ = 0;
    return *this;
  }

template<class T>
inline complex<T>& complex<T>::operator=(complex<float> const& rhs)
  {
    re_ = rhs.real();
    im_ = rhs.imag();
    return *this;
  }

template<class T>
inline complex<T>& complex<T>::operator=(complex<double> const& rhs)
  {
    re_ = rhs.real();
    im_ = rhs.imag();
    return *this;
  }

template<class T>
inline complex<T>& complex<T>::operator+=(T const& rhs)
  {
    re_ += rhs;
    return *this;
  }

template<class T>
inline complex<T>& complex<T>::operator+=(complex<float> const& rhs)
  {
    re_ += rhs.real();
    im_ += rhs.imag();
    return *this;
  }

template<class T>
inline complex<T>& complex<T>::operator+=(complex<double> const& rhs)
  {
    re_ += rhs.real();
    im_ += rhs.imag();
    return *this;
  }

template<class T>
inline complex<T>& complex<T>::operator-=(T const& rhs)
  {
    re_ -= rhs;
    return *this;
  }

template<class T>
inline complex<T>& complex<T>::operator-=(complex<float> const& rhs)
  {
    re_ -= rhs.real();
    im_ -= rhs.imag();
    return *this;
  }

template<class T>
inline complex<T>& complex<T>::operator-=(complex<double> const& rhs)
  {
    re_ -= rhs.real();
    im_ -= rhs.imag();
    return *this;
  }

template<class T>
inline complex<T>& complex<T>::operator*=(T const& rhs)
  {
    re_ *= rhs;
    im_ *= rhs;
    return *this;
  }

template<class T>
inline complex<T>& complex<T>::operator/=(T const& rhs)
  {
    re_ /= rhs;
    im_ /= rhs;
    return *this;
  }


// Implementation of complex<T> free fns

template<class T>
inline complex<T> operator+(complex<T> const& x, complex<T> const& y)
{
  complex<T> copy(x);
  return copy += y;
}

template<class T>
inline complex<T> operator+(complex<T> const& x, T const& y)
{
  complex<T> copy(x);
  return copy += y;
}

template<class T>
inline complex<T> operator+(T const& x, complex<T> const& y)
{
  complex<T> copy(y);
  return copy += x;
}

template<class T>
inline complex<T> operator-(complex<T> const& x, complex<T> const& y)
{
  complex<T> copy(x);
  return copy -= y;
}

template<class T>
inline complex<T> operator-(complex<T> const& x, T const& y)
{
  complex<T> copy(x);
  return copy -= y;
}

template<class T>
inline complex<T> operator-(T const& x, complex<T> const& y)
{
  complex<T> copy(-y);
  return copy += x;
}

template<class T>
inline complex<T> operator*(complex<T> const& x, complex<T> const& y)
{
  complex<T> copy(x);
  return copy *= y;
}

template<class T>
inline complex<T> operator*(complex<T> const& x, T const& y)
{
  complex<T> copy(x);
  return copy *= y;
}

template<class T>
inline complex<T> operator*(T const& x, complex<T> const& y)
{
  complex<T> copy(y);
  return copy *= x;
}

template<class T>
inline complex<T> operator/(complex<T> const& x, complex<T> const& y)
{
  complex<T> copy(x);
  return copy /= y;
}

template<class T>
inline complex<T> operator/(complex<T> const& x, T const& y)
{
  complex<T> copy(x);
  return copy /= y;
}

template<class T>
inline complex<T> operator+(complex<T> const& x)
{
  return x;
}

template<class T>
inline complex<T> operator-(complex<T> const& x)
{
  return complex<T>(-x.real(), -x.imag());
}

template<class T>
inline bool operator==(complex<T> const& x, complex<T> const& y)
{
  return x.real() == y.real() && x.imag() == y.imag();
}

template<class T>
inline bool operator==(complex<T> const& x, T const& y)
{
  return x.imag() == 0 && x.real() == y;
}

template<class T>
inline bool operator==(T const& x, complex<T> const& y)
{
  return y.imag() == 0 && x == y.real();
}

template<class T>
inline bool operator!=(complex<T> const& x, complex<T> const& y)
{
  return !(x == y);
}

template<class T>
inline bool operator!=(complex<T> const& x, T const& y)
{
  return !(x == y);
}

template<class T>
inline bool operator!=(T const& x, complex<T> const& y)
{
  return !(x == y);
}

template<class T>
inline T real(complex<T> const& x)
{
  return x.real();
}

template<class T>
inline T imag(complex<T> const& x)
{
  return x.imag();
}

template<class T>
inline T abs(complex<T> const& x)
{
  return sqrt(norm(x));
}

template<class T>
inline T arg(complex<T> const& x)
{
  return atan2(x.imag(), x.real());
}

template<class T>
inline T norm(complex<T> const& x)
{
  T re = x.real();
  T im = x.imag();
  return (re*re)+(im*im);
}

template<class T>
inline complex<T> conj(complex<T> const& x)
{
  return complex<T>(x.real(), -x.imag());
}

template<class T>
inline complex<T> polar(T const& mag, T const& theta)
{
  return complex<T>(mag*cos(theta), mag*sin(theta));
}

template<class T>
inline complex<T> cos(complex<T> const& x)
{
  return complex<T>(cos(x.real())*cosh(x.imag()), -sin(x.real())*sinh(x.imag()));
}

template<class T>
inline complex<T> cosh(complex<T> const& x)
{
  return complex<T>(cosh(x.real())*cos(x.imag()), sinh(x.real())*sin(x.imag()));
}

template<class T>
inline complex<T> exp(complex<T> const& x)
{
  return polar(T(::exp(x.real())), x.imag());
}

template<class T>
inline complex<T> log(complex<T> const& x)
{
  return complex<T>(log(abs(x)), arg(x));
}

template<class T>
inline complex<T> log10(complex<T> const& x)
{
  return log(x)/=log(10);
}

template<class T>
inline complex<T> pow(complex<T> const& x, T const& y)
{
  return exp(y*log(x));
}

template<class T>
inline complex<T> pow(complex<T> const& x, complex<T> const& y)
{
  return exp(y*log(x));
}

template<class T>
inline complex<T> pow(T const& x, complex<T> const& y)
{
  return exp(y*T(log(x)));
}

template<class T>
inline complex<T> sin(complex<T> const& x)
{
  return complex<T>(sin(x.real())*cosh (x.imag()), cos(x.real())*sinh (x.imag()));
}

template<class T>
inline complex<T> sinh(complex<T> const& x)
{
  return complex<T>(sinh(x.real())*cos(x.imag()), cosh(x.real())*sin(x.imag()));
}

template<class T>
inline complex<T> tan(complex<T> const& x)
{
  return sin(x) /= cos(x);
}

template<class T>
inline complex<T> tanh(complex<T> const& x)
{
  return sinh(x) /= cosh(x);
}

#endif
